{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 第19章 马尔可夫链蒙特卡罗法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 习题19.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;用蒙特卡罗积分法求\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} x^2 \\exp \\left( - \\frac{x^2}{2}\\right) \\mathrm{d} x\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "\n",
    "1. 给出蒙特卡罗积分计算方法\n",
    "2. 将被积分的函数分解，得到其数学期望表示形式\n",
    "3. 自编程实现，使用样本均值求近似计算积分"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：蒙特卡罗积分计算方法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19.1.3节的关于蒙特卡罗积分法的描述：\n",
    "\n",
    "> &emsp;&emsp;假设有一个函数 $h(x)$，目标是计算该函数的积分\n",
    "> \n",
    "> $$ \\int_{\\mathcal{X}} h(x) \\mathrm{d} x $$\n",
    "> \n",
    "> &emsp;&emsp;如果能够将函数 $h(x)$ 分解成一个函数 $f(x)$ 和一个概率密度函数 $p(x)$ 的乘积的形式, 那么就有\n",
    "> \n",
    "> $$ \\int_{\\mathcal{X}} h(x) \\mathrm{d} x = \\int_{\\mathcal{X}} f(x) p(x) \\mathrm{d} x = E_{p(x)}[f(x)] $$\n",
    "> \n",
    "> 于是函数 $h(x)$ 的积分可以表示为函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望。实际上, 给定一个概率密度函数 $p(x)$, 只要取 $\\displaystyle f(x)=\\frac{h(x)}{p(x)}$，就可得上式。就是说, 任何一个函数的积分都可以表示为某一个函数的数学期望的形式。而函数的数学期望又可以通过函数的样本均值估计。于是，就可以利用样本均值来近似计算积分。\n",
    "> \n",
    "> $$ \\int_{\\mathcal{X}} h(x) \\mathrm{d} x = E_{p(x)}[f(x)] \\approx \\frac{1}{n} \\sum_{i=1}^{n} f (x_i) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：将所积分的函数分解，其得到数学期望表示形式**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;取概率密度函数为标准正态分布的密度函数 $\\displaystyle p(x) = \\frac{1}{\\sqrt{2 \\pi}} \\exp \\left( - \\frac{x^2}{2} \\right)$，则\n",
    "$$\n",
    "f(x) = \\frac{h(x)}{p(x)} = \\frac{\\displaystyle x^2 \\exp \\left( - \\frac{x^2}{2} \\right)}{\\displaystyle \\frac{1}{\\sqrt{2 \\pi}} \\exp \\left(- \\frac{x^2}{2} \\right)} = \\sqrt{2 \\pi} x^{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;将原函数积分表示为函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望\n",
    "\n",
    "$$\n",
    "\\int_{\\mathcal{X}} x^2 \\exp \\left( - \\frac{x^2}{2} \\right) \\mathrm{d} x = \\int_{\\mathcal{X}} f(x) p(x) \\mathrm{d} x = E_{p(x)}[f(x)]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：自编程实现，使用样本均值求近似计算积分**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class MonteCarloIntegration:\n",
    "    def __init__(self, func_f, func_p):\n",
    "        # 所求期望的函数\n",
    "        self.func_f = func_f\n",
    "        # 抽样分布的概率密度函数\n",
    "        self.func_p = func_p\n",
    "\n",
    "    def solve(self, num_samples):\n",
    "        \"\"\"\n",
    "        蒙特卡罗积分法\n",
    "        :param num_samples: 抽样样本数量\n",
    "        :return: 样本的函数均值\n",
    "        \"\"\"\n",
    "        samples = self.func_p(num_samples)\n",
    "        vfunc_f = lambda x: self.func_f(x)\n",
    "        vfunc_f = np.vectorize(vfunc_f)\n",
    "        y = vfunc_f(samples)\n",
    "        return np.sum(y) / num_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "抽样样本数量: 1000000.0\n",
      "近似解: 2.5084502655966445\n"
     ]
    }
   ],
   "source": [
    "def func_f(x):\n",
    "    \"\"\"定义函数f\"\"\"\n",
    "    return x ** 2 * np.sqrt(2 * np.pi)\n",
    "\n",
    "\n",
    "def func_p(n):\n",
    "    \"\"\"定义在分布上随机抽样的函数g\"\"\"\n",
    "    return np.random.standard_normal(int(n))\n",
    "\n",
    "\n",
    "# 设置样本数量\n",
    "num_samples = 1e6\n",
    "\n",
    "# 使用蒙特卡罗积分法进行求解\n",
    "monte_carlo_integration = MonteCarloIntegration(func_f, func_p)\n",
    "result = monte_carlo_integration.solve(num_samples)\n",
    "print(\"抽样样本数量:\", num_samples)\n",
    "print(\"近似解:\", result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;证明如果马尔可夫链是不可约的，且有一个状态是非周期的，则其他所有状态也是非周期的，即这个马尔可夫链是非周期的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**  \n",
    "\n",
    "1. 给出马尔科夫链的基本定义\n",
    "2. 给出不可约的定义\n",
    "3. 给出非周期的定义\n",
    "4. 证明如果马尔可夫链是不可约的，则所有状态周期都相同\n",
    "5. 证明原命题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：马尔可夫链**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.1马尔科夫链的基本定义：\n",
    "\n",
    "> **定义19.1（马尔科夫链）** 考虑一个随机变量的序列$X = \\{ X_0, X_1, \\cdots, X_t, \\cdots \\}$，这里$X_t$表示时刻$t$的随机变量，$t=0,1,2,\\cdots$。每个随机变量$X_t(t=0, 1, 2, \\cdots)$的取值集合相同，称为状态空间，表示为$\\mathcal{S}$。随机变量可以是离散的，也可以是连续的。以上随机变量的序列构成随机过程（stochastic process）。  \n",
    "> \n",
    "> &emsp;&emsp;假设在时刻0的随机变量$X_0$遵循概率分布$P(X_0) = \\pi_0$，称为初始状态分布。在某个时刻$t \\geqslant 1$的随机变量$X_t$与前一时刻的随机变量$X_{t-1}$之间有条件分布$P(X_t | X_{t-1})$，如果$X_t$只依赖于$X_{t-1}$，而不依赖于过去的随机变量$\\{ X_0, X_1, \\cdots, X_{t-2}\\}$，这一性质称为马尔可夫性，即\n",
    "> \n",
    "> $$ P(X_t | X_0, X_1, \\cdots, X_{t-1}) = P(X_t | X_{t-1}), t=1,2,\\cdots $$\n",
    "> \n",
    "> 具有马尔可夫性的随机序列$X = \\{ X_0, X_1, \\cdots, X_t, \\cdots \\}$称为马尔可夫链，或马尔可夫过程。条件概率分布$P(X_t | X_{t-1})$称为马尔可夫链的转移概率分布。转移概率分布决定了马尔可夫链的特性。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19.2.2节的状态转移概率：\n",
    "\n",
    "> $$ P_{ij}^t = P(X_t = i | X_0 = j) $$\n",
    "> 表示时刻0从状态$j$出发，时刻$t$到达状态$i$的$t$步转移概率。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：不可约的定义**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.3的不可约的定义：\n",
    "\n",
    "> &emsp;&emsp;**定义19.3（不可约）**  设有马尔可夫链 $X= \\{ X_0, X_1, \\cdots, X_t, \\cdots \\}$，状态空间为 $\\mathcal{S}$，对于任意状态 $i, j \\in \\mathcal{S}$，如果存在一个时刻 $t(t>0)$ 满足\n",
    "> \n",
    "> $$ P (X_t = i | X_0 =j ) > 0 $$\n",
    "> \n",
    "> 也就是说，时刻0从状态 $j$ 出发, 时刻 $t$ 到达状态 $i$ 的概率大于0，则称此马尔可夫链 $X$ 是不可约的（irreducible），否则称马尔可夫链是可约的（reducible）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：非周期的定义**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.4的非周期的定义：\n",
    "\n",
    "> &emsp;&emsp;**定义19.4（非周期）** 设有马尔可夫链 $X = \\{X_0, X_1, \\cdots, X_t, \\cdots \\}$，状态空间为 $\\mathcal{S}$，对于任意状态 $i \\in \\mathcal{S}$，如果时刻0从状态 $i$ 出发，$t$ 时刻返回状态的所有时间长 $\\{ t: P( X_t = i | X_0 = i ) > 0 \\}$ 的最大公约数是1，则称此马尔可夫链 $X$ 是非周期的（aperiodic），否则称马尔可夫链是周期的（periodic）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第4步：证明如果马尔可夫链是不可约的，则所有状态周期都相同**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;如果该马尔可夫链是不可约的，则对于 $\\forall i, j \\in S$，必然存在时刻 $m, n (m > 0, n > 0)$ 满足 \n",
    "\n",
    "$$\n",
    "p_{ij}^{m} > 0, p_{ji}^{n}>0\n",
    "$$\n",
    "\n",
    "根据状态转移概率的定义，有 $p_{ij}^{m} = P ( X_m = i | X_0 = j), p_{ji}^{n} = P ( X_n = j | X_0 = i)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;可得\n",
    "\n",
    "$$\n",
    "p_{i i}^{m+n} = \\sum_{k \\in S} p_{i k}^{m} p_{k i}^{n} \\geqslant p_{i j}^{m} p_{j i}^{n}> 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;则必然存在 $s$ 使得 $P_{j j}^{s}>0$，且满足 \n",
    "\n",
    "$$\n",
    "P_{i i}^{m+n+s} \\geqslant P_{i j}^{m} P_{j j}^{s} P_{j i}^{n} > 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;假设状态 $i, j$ 的周期分别是 $d(i), d(j)$，由上两式可知\n",
    "\n",
    "$$\n",
    "(m+n) \\bmod d(i) = 0 \\\\\n",
    "(m+n+s) \\bmod d(i) = 0\n",
    "$$\n",
    "\n",
    "&emsp;&emsp;所以 $s \\bmod d(i) = 0$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;又因为 $P_{j j}^{s} > 0$，所以 $s \\bmod d(j) = 0$。因此 $d(i) \\bmod d(j) = 0$。同理可得 $d(j) \\bmod d(i) = 0$  \n",
    "\n",
    "&emsp;&emsp;故 $d(i)==d(j)$，可得如果马尔可夫链是不可约的，所有状态周期都是相同的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第5步：证明原命题**\n",
    "\n",
    "&emsp;&emsp;在不可约的马尔可夫链中如果有一个状态的是非周期的，即其周期为1，则其他所有状态的周期也为1，所以该马尔可夫链是非周期的。原命题得证。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;验证具有以下转移概率矩阵的马尔可夫链是可约的，但是非周期的。\n",
    "\n",
    "$$\n",
    "P=\\left[\\begin{array}{cccc}\n",
    "1 / 2 & 1 / 2 & 0 & 0 \\\\\n",
    "1 / 2 & 0 & 1 / 2 & 0 \\\\\n",
    "0 & 1 / 2 & 0 & 0 \\\\\n",
    "0 & 0 & 1 / 2 & 1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "\n",
    "1. 给出平稳分布的定义\n",
    "2. 计算该马尔可夫链的平稳分布\n",
    "3. 验证该马尔可夫链是可约的\n",
    "4. 验证该马尔可夫链是非周期的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：平稳分布的定义**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.2的平稳分布的定义：\n",
    "\n",
    "> &emsp;&emsp;**定义19.2（平稳分布）** 设有马尔可夫链 $X=\\{ X_0, X_1, \\cdots, X_t, \\cdots \\}$，其状态空间为 $\\mathcal{S}$，转移概率矩阵为 $P = (p_{i j})$，如果存在状态空间 $\\mathcal{S}$ 上的一个分布\n",
    "> $$ \\pi = \\left[ \\begin{array}{c}\n",
    "\\pi_1 \\\\\n",
    "\\pi_2 \\\\\n",
    "\\vdots\n",
    "\\end{array} \\right] $$\n",
    "> 使得\n",
    "> $$ \\pi = P \\pi $$\n",
    "> 则称 $\\pi$ 为马尔可夫链 $X = \\{X_0, X_1, \\cdots, X_t, \\cdots \\}$ 的平稳分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的引理19.1：\n",
    "\n",
    "> **引理19.1** 给定一个马尔可夫链 $ X = \\{X_0, X_1, \\cdots, X_t, \\cdots \\}$，状态空间为 $\\mathcal{S}$，转移概率矩阵为 $P = (p_{i j})$，则分布 $\\pi = (\\pi_1, \\pi_2, \\cdots )^T$ 为 $X$ 的平稳分布的充分必要条件是 $ \\pi = (\\pi_1, \\pi_2, \\cdots )^T$ 是下列方程组的解：\n",
    "> $$ \\begin{array}{l}\n",
    "\\displaystyle x_i = \\sum_j p_{i j} x_j, \\quad i = 1,2, \\cdots \\\\\n",
    "x_i \\geqslant 0, \\quad i = 1, 2, \\cdots \\\\\n",
    "\\displaystyle \\sum_i x_i = 1\n",
    "\\end{array} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：计算该马尔可夫链的平稳分布**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;设平稳分布为 $\\pi = (x_1, x_2, x_3, x_4)^T$，则由引理19.1可得\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\displaystyle x_1 = \\frac{1}{2} x_1 + \\frac{1}{2} x_2 \\\\\n",
    "\\displaystyle x_2 = \\frac{1}{2} x_1 + \\frac{1}{2} x_3 \\\\\n",
    "\\displaystyle x_3 = \\frac{1}{2} x_2 \\\\\n",
    "\\displaystyle x_4 = \\frac{1}{2} x_3 + 1 x_4 \\\\\n",
    "x_1 + x_2 + x_3 + x_4 = 1 \\\\\n",
    "x_i \\geqslant 0, \\quad i = 1, 2, 3, 4\n",
    "\\end{array}\n",
    "$$\n",
    "解方程组，得到唯一的平稳分布\n",
    "$$\n",
    "\\pi = (0, 0, 0, 1)^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：验证该马尔可夫链是不可约的**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;此马尔可夫链，转移到状态4后，就在该状态上循环跳转，不能到达状态1、状态2和状态3，最后停留在状态4，故该马尔可夫链是可约的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第4步：验证该马尔可夫链是非周期的**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 对于状态1，它返回状态1的最小时间长为1，故其所有时间长的最大公约数是1；\n",
    "2. 对于状态2，它通过$2 \\rightarrow 1 \\rightarrow 2$返回状态2的时间长为2，通过$2 \\rightarrow 1 \\rightarrow 1 \\rightarrow 2$返回状态2的时间长为3，故其所有时间长的最大公约数是1\n",
    "3. 对于状态3，它通过$3 \\rightarrow 2 \\rightarrow 3$返回状态3的时间长为2，通过$3 \\rightarrow 2 \\rightarrow 1 \\rightarrow 1 \\rightarrow 2 \\rightarrow 3$返回状态3的时间长为5，故其所有时间长的最大公约数是1\n",
    "4. 对于状态4，它返回状态4的最小时间长为1，故其所有时间长的最大公约数是1\n",
    "\n",
    "综上，该马尔可夫链是非周期的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;验证具有以下转移概率矩阵的马尔可夫链是不可约的，但是周期性的。\n",
    "\n",
    "$$\n",
    "P=\\left[\\begin{array}{cccc}\n",
    "0 & 1 / 2 & 0 & 0 \\\\\n",
    "1 & 0 & 1 / 2 & 0 \\\\\n",
    "0 & 1 / 2 & 0 & 1 \\\\\n",
    "0 & 0 & 1 / 2 & 0\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "1. 计算该马尔可夫链的平稳分布\n",
    "2. 验证该马尔可夫链是不可约的\n",
    "3. 验证该马尔可夫链是周期性的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：计算该马尔可夫链的平稳分布**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.2的平稳分布的定义和引理19.1（同习题19.4），设平稳分布为 $\\pi=(x_{1}, x_{2}, x_{3}, x_{4})^{\\mathrm{T}}$，则由引理19.1可得：\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\displaystyle x_1 = \\frac{1}{2} x_2 \\\\\n",
    "\\displaystyle x_2 = 1 x_1 + \\frac{1}{2} x_3 \\\\\n",
    "\\displaystyle x_3 = \\frac{1}{2} x_2 + 1 x_4 \\\\\n",
    "\\displaystyle x_4 = \\frac{1}{2} x_3 \\\\\n",
    "x_1 + x_2 + x_3 + x_4 = 1 \\\\\n",
    "x_i \\geqslant 0, \\quad i = 1, 2, 3, 4\n",
    "\\end{array}\n",
    "$$\n",
    "解方程组，得到唯一的平稳分布\n",
    "$$\n",
    "\\pi=(\\frac{1}{6}, \\frac{1}{3}, \\frac{1}{3}, \\frac{1}{6})^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：验证该马尔可夫链是不可约的**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;观察到平稳分布各项均大于0，说明此马尔可夫链从任意状态出发，经过若干步之后，可以到达任意状态，故该马尔可夫链是不可约的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：验证该马尔可夫链是周期性的**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;该马尔可夫链的转移仅发生在相邻奇偶状态之间，从每个状态出发，返回该状态的时刻都是2的倍数：$\\{2, 4, 6, \\cdots, 2n\\}, n \\in N^*$。故该马尔可夫链是周期性的，其周期为2。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;证明可逆马尔可夫链一定是不可约的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路**\n",
    "\n",
    "1. 给出可逆马尔可夫链的定义\n",
    "2. 构造特殊的可逆马尔可夫链\n",
    "3. 验证其是可约的，原命题不成立"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：可逆马尔可夫链**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的定义19.6的可逆马尔可夫链的定义：\n",
    "\n",
    "> &emsp;&emsp;**定义19.6（可逆马尔可夫链）** 设有马尔可夫链 $ X = \\{X_0, X_1, \\cdots, X_t, \\cdots \\}$，状态空间为 $\\mathcal{S}$，转移概率矩阵为 $P$，如果有状态分布 $\\pi = (\\pi_1, \\pi_2, \\cdots )^T$，对于任意状态 $i, j \\in \\mathcal{S}$，对任意一个时刻 $t$ 满足\n",
    "> $$ P (X_t = i | X_{t-1} = j) \\pi_j = P (X_{t-1} = j | X_t = i) \\pi_i, \\quad i,j = 1,2, \\cdots $$\n",
    "> 或简写为\n",
    "> $$ p_{j i} \\pi_j = p_{i j} \\pi_{i}, \\quad i, j=1,2, \\cdots $$\n",
    "> 则称此马尔可夫链 $X$ 为可逆马尔可夫链，上式称为细致平衡方程。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：构造特殊的可逆马尔可夫链**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;对于以单位矩阵做转移矩阵的马尔可夫链，其转移矩阵P\n",
    "$$\n",
    "P = \\left [ \\begin{array}{cccc}\n",
    "1 & 0 & \\cdots & 0 \\\\\n",
    "0 & 1 & \\cdots & 0 \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "0 & 0 & 0 & 1\n",
    "\\end{array} \\right ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;对于任意平稳分布 $\\pi$，满足$P \\pi = \\pi$ 恒成立，所以该马尔可夫链是可逆的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：验证其是可约的，原命题被证伪**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;对于任意状态 $i \\in S$，跳转到其他的状态 $j \\in S(j \\neq i)$ 的概率始终为0，即 $p^t_{i, j} = 0, \\forall t \\in N^{*}$。所以该马尔可夫链是可约的。\n",
    "\n",
    "&emsp;&emsp;这与原命题矛盾，故原命题不成立。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;从一般的Metropolis-Hastings算法推导出单分量Metropolis-Hastings算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "\n",
    "1. 给出一般的Metropolis-Hastings算法\n",
    "2. 给出单分量Metropolis-Hastings算法的基本思想\n",
    "3. 推导单分量Metropolis-Hastings算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：一般的Metropolis-Hastings算法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19章的算法19.2的Metropolis-Hastings算法：\n",
    "\n",
    "> **算法19.2（Metropolis-Hastings算法）**  \n",
    "> 输入：抽样的目标分布的密度函数$p(x)$，函数$f(x)$；  \n",
    "> 输出：$p(x)$的随机样本$x_{m+1}, x_{m+2}, \\cdots, x_n$，函数样本均值$f_{m n}$；  \n",
    "> 参数：收敛步数$m$，迭代步数$n$。  \n",
    "> （1）任意选择一个初始值$x_0$  \n",
    "> （2） 对 $i=1,2, \\cdots, n$循环执行  \n",
    "> &emsp;&emsp;（a）设状态$x_{i-1} = x$，按照建议分布$q (x, x')$随机抽取一个候选状态$x'$。  \n",
    "> &emsp;&emsp;（b）计算接受概率\n",
    "> $$ \\alpha (x, x') = \\min \\left \\{ 1, \\frac{p (x') q (x', x)}{p(x) q(x, x')} \\right \\} $$\n",
    "> &emsp;&emsp;（c）从区间$(0,1)$中按均匀分布随机抽取一个数$u$。  \n",
    "> &emsp;&emsp;&emsp;&emsp;&nbsp;若 $u \\leqslant \\alpha (x, x')$，则状态 $x_i = x'$；否则，状态 $x_i = x$。  \n",
    "> （3）得到样本集合 $ \\{ x_{m+1}, x_{m+2}, \\cdots, x_n \\}$  \n",
    "> 计算\n",
    "> $$ f_{m n} = \\frac{1}{n-m} \\sum_{i = m + 1}^n f(x_i) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. 理解单分量 Metropolis-Hastings 算法的基本思想**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19.4.3节的单分量Metropolis-Hastings算法：\n",
    "\n",
    "> &emsp;&emsp;在 Metropolis-Hastings 算法中，通常需要对多元变量分布进行抽样，有时对多元变量分布的抽样是困难的。可以对多元变量的每一变量的条件分布依次分别进行抽样，从而实现对整个多元变量的一次抽样, 这就是单分量 Metropolis-Hastings 算法。  \n",
    "> &emsp;&emsp;假设马尔可夫链的状态由 $k$ 维随机变量表示\n",
    "> $$ x = (x_1, x_2, \\cdots, x_k)^T $$\n",
    "> 其中 $x_{j}$ 表示随机变量 $x$ 的第 $j$ 个分量，$j=1,2, \\cdots, k$，而 $x^{(i)}$ 表示马尔可夫链在时刻 $i$ 的状态\n",
    "> $$ x^{(i)} = (x_1^{(i)}, x_2^{(i)}, \\cdots, x_k^{(i)})^T, \\quad i = 1, 2, \\cdots, n $$\n",
    "> 其中 $x_j^{(i)}$ 是随机变量 $x^{(i)}$ 的第 $j$ 个分量，$j = 1, 2, \\cdots, k$。  \n",
    "> &emsp;&emsp;为了生成容量为 $n$ 的样本集合 $\\{ x^{(1)}, x^{(2)}, \\cdots, x^{(n)} \\}$，单分量 Metropolis-Hastings 算法由下面的 $k$ 步迭代实现 Metropolis-Hastings 算法的一次迭代。  \n",
    "> &emsp;&emsp;设在第 $(i-1)$ 次迭代结束时分量 $x_j$ 的取值为 $x_j^{(i-1)}$，在第 $i$ 次迭代的第 $j$ 步，对分量 $x_j$ 根据 Metropolis-Hastings 算法更新，得到其新的取值 $x_j^{(i)}$。首先，由建议分布 $q (x_j^{(i-1)}, x_j | x_{-j}^{(i)})$ 抽样产生分量 $x_j$ 的候选值 $x{'}_j^{(i)}$，这里 $x_{-j}^{(i)}$ 表示在第 $i$ 次迭代的第 $(j-1)$ 步后的 $x^{(i)}$ 除去 $x_j^{(i-1)}$ 的所有值，即\n",
    "> $$ x_{-j}^{(i)}=(x_1^{(i)}, \\cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \\cdots, x_k^{(i-1)})^T $$\n",
    "> 其中分量 $1, 2, \\cdots, j-1$ 已经更新。然后，按照接受概率\n",
    "> $$ \\alpha \\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right ) = \\min \\left\\{ 1, \\frac\n",
    "{p \\left( x{'}_j^{(i)} | x_{-j}^{(i)} \\right) q \\left( x{'}_j^{(i)}, x_j^{(i-1)} | x_{-j}^{(i)} \\right)}\n",
    "{p \\left( x_j^{(i-1)} | x_{-j}^{(i)} \\right) q \\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right)} \n",
    "\\right \\} $$\n",
    "> 抽样决定是否接受候选值 $x{'}_j^{(i)}$ 。如果 $x{'}_j^{(i)}$ 被接受，则令 $x_j^{(i)}=x{'}_j^{(i)}$；否则令 $x_j^{(i)} = x_{j}^{(i-1)}$。其余分量在第 $j$ 步不改变。马尔可夫链的转移概率为\n",
    "> $$ p\\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right) \n",
    "= \\alpha \\left(\n",
    "x_j^{(i-1)}, x{'}_{j}^{(i)} | x_{-j}^{(i)} \\right) q \\left(x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3. 推导单分量 Metropolis-Hastings 算法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**单分量 Metropolis-Hastings 算法**  \n",
    "输入：抽样的目标分布的密度函数 $p(x)$，函数 $f(x)$，其中 $x = \\left( x_1, x_2, \\cdots, x_k \\right)^T$；  \n",
    "输出：$p(x)$ 的随机样本 $\\{ x^{(m+1)}, x^{(m+2)}, \\cdots, x^{(n)} \\}$，函数样本均值 $f_{m n}$；  \n",
    "参数：收敛步数 $m$, 迭代步数 $n$。  \n",
    "（1）任意选择一个初始值 $x_0$  \n",
    "（2）对 $i = 1, 2, \\cdots, n$ 循环执行  \n",
    "&emsp;&emsp;设在第 $(i-1)$ 次迭代结束时分量 $x_{j}$ 的取值为 $x_{j}^{(i-1)}$  \n",
    "&emsp;&emsp;（a）对 $j=1,2, \\cdots, k$ 循环执行  \n",
    "&emsp;&emsp;&emsp;&emsp;（i）在第 $i$ 次迭代的第 $j$ 步，由建议分布 $q (x_j^{(i-1)}, x_j | x_{-j}^{(i)})$ 抽样产生分量 $x_j$ 的候选值 $x{'}_j^{(i)}$，  \n",
    "&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;&nbsp;&nbsp;这里 $x_{-j}^{(i)}$ 表示在第 $i$ 次迭代的第 $(j-1)$ 步后的 $x^{(i)}$ 除去 $x_j^{(i-1)}$ 的所有值，即\n",
    "$$\n",
    "x_{-j}^{(i)}=(x_1^{(i)}, \\cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \\cdots, x_k^{(i-1)})^T\n",
    "$$  \n",
    "&emsp;&emsp;&emsp;&emsp;（ii）计算接受概率\n",
    "$$\n",
    "\\alpha \\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right ) = \\min \\left\\{ 1, \\frac\n",
    "{p \\left( x{'}_j^{(i)} | x_{-j}^{(i)} \\right) q \\left( x{'}_j^{(i)}, x_j^{(i-1)} | x_{-j}^{(i)} \\right)}\n",
    "{p \\left( x_j^{(i-1)} | x_{-j}^{(i)} \\right) q \\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right)} \n",
    "\\right \\}\n",
    "$$\n",
    "&emsp;&emsp;&emsp;&emsp;（iii）从区间 $(0,1)$ 中按均匀分布随机抽取一个数 $u$。  \n",
    "&emsp;&emsp;&emsp;&emsp;&emsp;&emsp;&nbsp;&nbsp;若 $u \\leqslant \\alpha \\left( x_j^{(i-1)}, x{'}_j^{(i)} | x_{-j}^{(i)} \\right )$，则状态 $x_j^{(i)} = x{'}_j^{(i)}$；否则，状态 $x_j^{(i)} = x_j^{(i-1)}$。  \n",
    "（3）得到样本集合 $\\{ x^{(m+1)}, x^{(m+2)}, \\cdots, x^{(n)} \\}$  \n",
    "计算\n",
    "$$\n",
    "f_{m n} = \\frac{1}{n - m} \\sum_{i = m + 1}^n f (x^{(i)})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;假设进行伯努利实验，后验概率为 $P(\\theta | y)$，其中变量 $y \\in \\{ 0, 1 \\}$ 表示实验可能的结果，变量 $\\theta$ 表示结果为 1 的概率。再假设先验概率 $P(\\theta)$ 遵循 Beta 分布 $B(\\alpha, \\beta)$，其中 $\\alpha = 1, \\beta = 1$；似然函数 $P(y | \\theta)$ 遵循二项分布$ \\text{Bin}(n, k, \\theta)$，其中 $n = 10, k = 4$，即实验进行 10 次其中结果为 1 的次数为 4。试用 Metropolis-Hastings 算法求后验概率分布 $P(\\theta | y) \\propto P(\\theta) P(y | \\theta)$ 的均值和方差。（提示：可采用 Metropolis 选择，即假设建议分布是对称的。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解题思路**\n",
    "\n",
    "1. 给出 Metropolis-Hastings 算法\n",
    "2. 写出目标分布和建议分布\n",
    "3. 自编程实现 Metropolis-Hastings 算法求解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：Metropolis-Hastings 算法**  \n",
    "\n",
    "&emsp;&emsp;Metropolis-Hastings 算法的步骤详见书中第373~374页算法19.2。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：写出目标分布和建议分布**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据题意，可知后验概率 $P(\\theta) \\propto B(1, 1)$，后验函数 $P(y | \\theta) \\propto \\text{Bin}(10, 4, \\theta)$\n",
    "\n",
    "&emsp;&emsp;则后验概率分布\n",
    "$$\n",
    "\\begin{aligned}\n",
    "P(\\theta | y)\n",
    "& \\propto P(\\theta) P(y | \\theta) \\\\\n",
    "& \\propto B(1, 1) \\text{Bin}(10, 4, \\theta)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "&emsp;&emsp;可取建议分布为 $B(1, 1)$，接受分布为 $\\text{Bin}(10, 4, \\theta)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：自编程实现 Metropolis-Hastings 算法求解**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import beta, binom\n",
    "\n",
    "\n",
    "class MetropolisHastings:\n",
    "    def __init__(self, proposal_dist, accepted_dist, m=1e4, n=1e5):\n",
    "        \"\"\"\n",
    "        Metropolis Hastings\n",
    "\n",
    "        :param proposal_dist: 建议分布\n",
    "        :param accepted_dist: 接受分布\n",
    "        :param m: 收敛步数\n",
    "        :param n: 迭代步数\n",
    "        \"\"\"\n",
    "        self.proposal_dist = proposal_dist\n",
    "        self.accepted_dist = accepted_dist\n",
    "        self.m = m\n",
    "        self.n = n\n",
    "\n",
    "    @staticmethod\n",
    "    def __calc_acceptance_ratio(q, p, x, x_prime):\n",
    "        \"\"\"\n",
    "        计算接受概率\n",
    "\n",
    "        :param q: 建议分布\n",
    "        :param p: 接受分布\n",
    "        :param x: 上一状态\n",
    "        :param x_prime: 候选状态\n",
    "        \"\"\"\n",
    "        prob_1 = p.prob(x_prime) * q.joint_prob(x_prime, x)\n",
    "        prob_2 = p.prob(x) * q.joint_prob(x, x_prime)\n",
    "        alpha = np.min((1., prob_1 / prob_2))\n",
    "        return alpha\n",
    "\n",
    "    def solve(self):\n",
    "        \"\"\"\n",
    "        Metropolis Hastings 算法求解\n",
    "        \"\"\"\n",
    "        all_samples = np.zeros(self.n)\n",
    "        # (1) 任意选择一个初始值\n",
    "        x_0 = np.random.random()\n",
    "        for i in range(int(self.n)):\n",
    "            x = x_0 if i == 0 else all_samples[i - 1]\n",
    "            # (2.a) 从建议分布中抽样选取\n",
    "            x_prime = self.proposal_dist.sample()\n",
    "            # (2.b) 计算接受概率\n",
    "            alpha = self.__calc_acceptance_ratio(self.proposal_dist, self.accepted_dist, x, x_prime)\n",
    "            # (2.c) 从区间 (0,1) 中按均匀分布随机抽取一个数 u\n",
    "            u = np.random.uniform(0, 1)\n",
    "            # 根据 u <= alpha，选择 x 或 x_prime 进行赋值\n",
    "            if u <= alpha:\n",
    "                all_samples[i] = x_prime\n",
    "            else:\n",
    "                all_samples[i] = x\n",
    "\n",
    "        # (3) 随机样本集合\n",
    "        samples = all_samples[self.m:]\n",
    "        # 函数样本均值\n",
    "        dist_mean = samples.mean()\n",
    "        # 函数样本方差\n",
    "        dist_var = samples.var()\n",
    "        return samples[self.m:], dist_mean, dist_var\n",
    "\n",
    "    @staticmethod\n",
    "    def visualize(samples, bins=50):\n",
    "        \"\"\"\n",
    "        可视化展示\n",
    "        :param samples: 抽取的随机样本集合\n",
    "        :param bins: 频率直方图的分组个数\n",
    "        \"\"\"\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.set_title('Metropolis Hastings')\n",
    "        ax.hist(samples, bins, alpha=0.7, label='Samples Distribution')\n",
    "        ax.set_xlim(0, 1)\n",
    "        ax.legend()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProposalDistribution:\n",
    "    \"\"\"\n",
    "    建议分布\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod\n",
    "    def sample():\n",
    "        \"\"\"\n",
    "        从建议分布中抽取一个样本\n",
    "        \"\"\"\n",
    "        # B(1,1)\n",
    "        return beta.rvs(1, 1, size=1)\n",
    "\n",
    "    @staticmethod\n",
    "    def prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        return beta.pdf(x, 1, 1)\n",
    "\n",
    "    def joint_prob(self, x_1, x_2):\n",
    "        \"\"\"\n",
    "        P(X = x_1, Y = x_2) 的联合概率\n",
    "        \"\"\"\n",
    "        return self.prob(x_1) * self.prob(x_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class AcceptedDistribution:\n",
    "    \"\"\"\n",
    "    接受分布\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod\n",
    "    def prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        # Bin(4, 10)\n",
    "        return binom.pmf(4, 10, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "均值: 0.4159270490443074\n",
      "方差: 0.019403228732637397\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfq0lEQVR4nO3de5RU5Z3u8e8jKERFRGgNt9iojYKgxDRq0BAMSUTjiMdIxBkVDZEYUJOsREVPljqOzJDRZaJHOaZP9HQblUuMF8Y5k2jwFuMFARkFFUVB6ECgQRERQcHf+aN2QwHddHV1VXfT+/ms1atqv/vdu361xad2vbUvigjMzCxd9mrpAszMrPk5/M3MUsjhb2aWQg5/M7MUcvibmaWQw9/MLIUc/pZakiol3ZQ8/5qkRS1dU0MkbZB0WEvXYXs+h7/lRdJSSZ9K6rZT+3xJIak0h3UMk1RdtCIbISL+EhFHNna57A+QrLbSZBu0b0pNkp6W9IOd6tw/It5tynrNwOFvTbMEOK92QtJA4AuFfIGmBqiZ1c3hb03xO+DCrOkxwL3ZHSR1kHSLpGWSVkm6S9IXJO0H/BfQIxnK2CCph6QbJD0o6T5J64GLkvaZkt6XtFjSJVnrr+0/XdJHkuZJOjZrfr9kD3qdpIWSzqzrjez8LUTS1ZL+lqxzkaTh+W4kSd+R9Iqk9ZKWS7oha17H5L2uTWp8WdIhkiYBXwPuSLbNHUn/kHRE8rxS0p2S/jOp8yVJh2et+9tJ7R9KmiLpmdpvEpKOSKY/lLRG0vR835/tmRz+1hQvAgckAdsOOBe4b6c+vwT6AoOAI4CewHUR8TFwGrAiGcrYPyJWJMuMBB4EDgTuB6YC1UAP4BzgX3cK45HA74GDgAeARyTtLWlv4D+Ax4GDgcuB+yXtdngnmX8ZMDgiOgGnAksbsV129jGZD8kDge8AP5J0VjJvDNAZ6A10BS4FPomI/wn8Bbgs2TaX1bPu84B/BroAi4FJyXvoRmYbXpOsdxEwJGu5fyGzXboAvYD/1YT3Z3sgh781Ve3e/7eAN4G/1c6QJOAS4KcR8X5EfAT8KzC6gXW+EBGPRMTnQDfgZODqiNgUEfOB3wIXZPWfGxEPRsRnwK1AR+DE5G9/YHJEfBoRTwKPkTVUVY+tQAegv6S9I2JpRLyzm/4/T/ba10laB7yaPTMino6I1yLi84h4lcyH2deT2Z+RCecjImJrRMyNiPUN1JftoYiYHRFbyHxQDkraTwcWRsRDybzbgb9nLfcZcCjQI9muzzXiNa0NcPhbU/0O+EfgInYa8gFKgH2BuVnB+MekfXeWZz3vAdR+cNR6j8w3iF36Jx8Ytd8SegDLk7b6lt1FRCwGfgLcAKyWNE1Sj90scktEHFj7BxyTPVPSCZKeklQj6UMye/e1P5T/DvgTME3SCkn/nnxjyVV2oG8k82EHyXvPek9BZrvUugoQMDsZDvt+I17T2gCHvzVJRLxH5off04GHdpq9BvgEODorHDtHRG1A1XdJ2ez2FcBBkjpltX2JrG8YZIZMAJC0F5lhjBXJX++krb5l63tfD0TEyWT2joPM8FW+HgBmAr0jojNwF5ngJSI+i4h/joj+ZIZlzmD77yhNueTuSjLbAdj2LWzbdET8PSIuiYgewA+BKbW/JVg6OPytEMYC30jG8bdJ9rj/D/ArSQcDSOop6dSkyyqgq6TO9a04IpYDzwP/lvw4ekzyevdndfuKpLOTI4N+Amwm83vES2TG269KfgMYBvwDMG13b0bSkZK+IakDsInMB9jWhjdDvTqR+faySdLxZL4p1b7WKZIGJr+ZrCczHFP7WquAfI/p/09goKSzku0yAfhi1uuOklT7YfABmQ+aprxH28M4/K3JIuKdiJhTz+yryfwQ+WJy9M6fgSOT5d4kM/79bjIsVN/QynlAKZk9+YeB6yPiiaz5j5L5sfkDMr8FnJ3sUX8KnEnmh+U1wBTgwuR1d6cDMDlZ5u9kfiy+toFldmc8cKOkj4DrgBlZ875I5ofZ9cAbwDNs/9H8NuAcSR9Iur0xLxgRa4BRwL8Da4H+wBwyH4wAg4GXJG0g863kxxGxJI/3Znso+WYutidLDps8IiLOb+laWrNk6Ksa+KeIeKql67GW5z1/szZK0qmSDkyGr64l8zvDiy1clrUSDn+ztuurwDtkhq/+ATgrIj5p2ZKstfCwj5lZCnnP38wshVrFRbO6desWpaWlLV2GmdkeZe7cuWsioqGTJuvUKsK/tLSUOXPqO1LQzMzqIum9fJf1sI+ZWQo5/M3MUsjhb2aWQq1izN8sjT777DOqq6vZtGlTS5dirVzHjh3p1asXe+/dmAu+7p7D36yFVFdX06lTJ0pLS8lcdNNsVxHB2rVrqa6upk+fPgVbr4d9zFrIpk2b6Nq1q4PfdksSXbt2Lfg3RIe/WQty8FsuivHvxOFvZpZCHvM3ayXGVr5c0PXdfdHgBvtMmjSJBx54gHbt2rHXXnvxm9/8hhNOOKGgdWQbNmwYt9xyC+Xl5Xmvo7KykiuvvJJevXqxYcMGDjvsMK6//nqGDMncn/66665j6NChfPOb36xz+UceeYS+ffvSv3//Ouffdddd7Lvvvlx44YWNrnfdunU88MADjB8/HoAVK1ZwxRVX8OCDD+bxTovL4W9Nlm9o5RJOVjwvvPACjz32GPPmzaNDhw6sWbOGTz/9tKXLysm5557LHXfcAcBTTz3F2WefzVNPPUW/fv248cYbd7vsI488whlnnFFn+G/ZsoVLL70077rWrVvHlClTtoV/jx49WmXwg4d9zFJr5cqVdOvWjQ4dOgDQrVs3evTI3EztxhtvZPDgwQwYMIBx48ZRe/XfYcOG8dOf/pShQ4fSr18/Xn75Zc4++2zKysr4xS9+AcDSpUs56qijGDNmDMcccwznnHMOGzdu3OX1H3/8cb761a9y3HHHMWrUKDZs2ADAxIkT6d+/P8cccww///nPG3wfp5xyCuPGjaOiogKAiy66aFvg7ryu559/npkzZ3LllVcyaNAg3nnnHYYNG8a1117L17/+dW677TZuuOEGbrnllm3rv++++xgyZAgDBgxg9uzZALv0GTBgAEuXLmXixIm88847DBo0iCuvvJKlS5cyYMAAIPMD/8UXX8zAgQP58pe/zFNPZe6pU1lZydlnn82IESMoKyvjqquuyvU/YZM4/M1S6tvf/jbLly+nb9++jB8/nmeeeWbbvMsuu4yXX36ZBQsW8Mknn/DYY49tm7fPPvvw7LPPcumllzJy5EjuvPNOFixYQGVlJWvXrgVg0aJFjBs3jldffZUDDjiAKVOm7PDaa9as4aabbuLPf/4z8+bNo7y8nFtvvZX333+fhx9+mIULF/Lqq69u+0BpyHHHHcebb+54d8661jVkyBDOPPNMbr75ZubPn8/hhx8OZPbYn3nmGX72s5/tsu6PP/6Y559/nilTpvD9739/t3VMnjyZww8/nPnz53PzzTfvMO/OO+8E4LXXXmPq1KmMGTNm2xE88+fPZ/r06bz22mtMnz6d5cuX5/S+m8Lhb5ZS+++/P3PnzqWiooKSkhLOPfdcKisrgcxQygknnMDAgQN58sknWbhw4bblzjzzTAAGDhzI0UcfTffu3enQoQOHHXbYttDq3bs3J510EgDnn38+zz333A6v/eKLL/L6669z0kknMWjQIKqqqnjvvfc44IAD6NixIz/4wQ946KGH2HfffXN6L3Xdl6Qx6zr33HPrnXfeeecBMHToUNavX8+6detyqmlnzz33HBdccAEARx11FIceeihvvfUWAMOHD6dz58507NiR/v378957eV+vLWce8zdLsXbt2jFs2DCGDRvGwIEDqaqqYvTo0YwfP545c+bQu3dvbrjhhh2OMa8dJtprr722Pa+d3rJlC7DroYk7T0cE3/rWt5g6deouNc2ePZtZs2Yxbdo07rjjDp588skG38crr7xCv379dmhr3759zuvab7/96l13Xe+lffv2fP7559vacjkGf3c3zsreju3atdu2HYvJe/5mKbVo0SLefvvtbdPz58/n0EMP3RZk3bp1Y8OGDXn9YLls2TJeeOEFAKZOncrJJ5+8w/wTTzyRv/71ryxevBiAjRs38tZbb7FhwwY+/PBDTj/9dH79618zf/78Bl/rmWeeoaKigksuuWSH9vrW1alTJz766KOc38v06dOBzJ57586d6dy5M6WlpcybNw+AefPmsWTJkgbXPXToUO6//34A3nrrLZYtW8aRRx6Zcx2F5j1/s1aiuY9+2rBhA5dffjnr1q2jffv2HHHEEVRUVHDggQdyySWXMHDgQEpLSxk8uPF19evXj6qqKn74wx9SVlbGj370ox3ml5SUUFlZyXnnncfmzZsBuOmmm+jUqRMjR45k06ZNRAS/+tWv6lz/9OnTee6559i4cSN9+vThD3/4wy57/h999FGd6xo9ejSXXHIJt99+e04fbF26dGHIkCGsX7+ee+65B4Dvfve73HvvvQwaNIjBgwfTt29fALp27cpJJ53EgAEDOO2005gwYcK29YwfP55LL72UgQMH0r59eyorK3fY429ureIevuXl5eGbuey5mvNQz3xeq7UeUvrGG2/sElhtwdKlSznjjDNYsGBBS5fSptT170XS3IjI66QJD/uYmaWQw9/MCqq0tNR7/XuAnMb8Jf0U+AEQwGvAxcC+wHSgFFgKfC8iPkj6XwOMBbYCV0TEnwpduO35Cn05gz1RRPjibtagYgzPN7jnL6kncAVQHhEDgHbAaGAiMCsiyoBZyTSS+ifzjwZGAFMktSt45WZ7uI4dO7J27dqi/I9tbUft9fw7duxY0PXmerRPe+ALkj4js8e/ArgGGJbMrwKeBq4GRgLTImIzsETSYuB44IXClW225+vVqxfV1dXU1NS0dCnWytXeyauQGgz/iPibpFuAZcAnwOMR8bikQyJiZdJnpaSDk0V6Ai9mraI6aduBpHHAOIAvfelLTXsXZnugvffeu6B3ZjJrjAbDX1IXMnvzfYB1wO8lnb+7Repo2+V7bURUABWQOdQzl2Kt+DwOb5YOuRzt801gSUTURMRnwEPAEGCVpO4AyePqpH810Dtr+V5khonMzKyVyCX8lwEnStpXmcMShgNvADOBMUmfMcCjyfOZwGhJHST1AcqA2YUt28zMmiKXMf+XJD0IzAO2AK+QGa7ZH5ghaSyZD4hRSf+FkmYAryf9J0TE1iLVb2ZmecjpaJ+IuB64fqfmzWS+BdTVfxIwqWmlmZlZsfgMXzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYp5PA3M0uhXO/ha7bHyufuZHdfNLgIlZi1Ht7zNzNLoQbDX9KRkuZn/a2X9BNJB0l6QtLbyWOXrGWukbRY0iJJpxb3LZiZWWM1GP4RsSgiBkXEIOArwEbgYWAiMCsiyoBZyTSS+gOjgaOBEcAUSe2KU76ZmeWjscM+w4F3IuI9YCRQlbRXAWclz0cC0yJic0QsARYDxxegVjMzK5DGhv9oYGry/JCIWAmQPB6ctPcElmctU5207UDSOElzJM2pqalpZBlmZtYUOYe/pH2AM4HfN9S1jrbYpSGiIiLKI6K8pKQk1zLMzKwAGrPnfxowLyJWJdOrJHUHSB5XJ+3VQO+s5XoBK5paqJmZFU5jwv88tg/5AMwExiTPxwCPZrWPltRBUh+gDJjd1ELNzKxwcjrJS9K+wLeAH2Y1TwZmSBoLLANGAUTEQkkzgNeBLcCEiNha0KrNzKxJcgr/iNgIdN2pbS2Zo3/q6j8JmNTk6szMrCh8hq+ZWQo5/M3MUsjhb2aWQr6qp1kd8rkSKPhqoLbn8J6/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGf5NWG5Xuikpm1fd7zNzNLIYe/mVkK5RT+kg6U9KCkNyW9Iemrkg6S9ISkt5PHLln9r5G0WNIiSacWr3wzM8tHrnv+twF/jIijgGOBN4CJwKyIKANmJdNI6g+MBo4GRgBTJLUrdOFmZpa/BsNf0gHAUOBugIj4NCLWASOBqqRbFXBW8nwkMC0iNkfEEmAxcHxhyzYzs6bIZc//MKAG+L+SXpH0W0n7AYdExEqA5PHgpH9PYHnW8tVJ2w4kjZM0R9KcmpqaJr0JMzNrnFzCvz1wHPC/I+LLwMckQzz1UB1tsUtDREVElEdEeUlJSU7FmplZYeQS/tVAdUS8lEw/SObDYJWk7gDJ4+qs/r2zlu8FrChMuWZmVggNhn9E/B1YLunIpGk48DowExiTtI0BHk2ezwRGS+ogqQ9QBswuaNVmZtYkuZ7hezlwv6R9gHeBi8l8cMyQNBZYBowCiIiFkmaQ+YDYAkyIiK0Fr9zMzPKWU/hHxHygvI5Zw+vpPwmYlH9ZZmZWTD7D18wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLoVwv7GZmORhb+XKjl7n7osFFqMRs97znb2aWQg5/M7MUcvibmaWQw9/MLIVyCn9JSyW9Jmm+pDlJ20GSnpD0dvLYJav/NZIWS1ok6dRiFW9mZvlpzJ7/KRExKCJq7+g1EZgVEWXArGQaSf2B0cDRwAhgiqR2BazZzMyaqCnDPiOBquR5FXBWVvu0iNgcEUuAxcDxTXgdMzMrsFzDP4DHJc2VNC5pOyQiVgIkjwcn7T2B5VnLVidtO5A0TtIcSXNqamryq97MzPKS60leJ0XECkkHA09IenM3fVVHW+zSEFEBVACUl5fvMt/MzIonpz3/iFiRPK4GHiYzjLNKUneA5HF10r0a6J21eC9gRaEKNjOzpmsw/CXtJ6lT7XPg28ACYCYwJuk2Bng0eT4TGC2pg6Q+QBkwu9CFm5lZ/nIZ9jkEeFhSbf8HIuKPkl4GZkgaCywDRgFExEJJM4DXgS3AhIjYWpTqzcwsLw2Gf0S8CxxbR/taYHg9y0wCJjW5OjMzKwqf4WtmlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczSyGHv5lZCvkG7nuAfG4Kbma2O97zNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlkMPfzCyFcg5/Se0kvSLpsWT6IElPSHo7eeyS1fcaSYslLZJ0ajEKNzOz/DXmJK8fA28AByTTE4FZETFZ0sRk+mpJ/YHRwNFAD+DPkvr6bl5mdcvnJL67LxpchEosTXLa85fUC/gO8Nus5pFAVfK8Cjgrq31aRGyOiCXAYjI3fDczs1Yi12GfXwNXAZ9ntR0SESsBkseDk/aewPKsftVJ2w4kjZM0R9KcmpqaxtZtZmZN0GD4SzoDWB0Rc3Ncp+poi10aIioiojwiyktKSnJctZmZFUIuY/4nAWdKOh3oCBwg6T5glaTuEbFSUndgddK/GuidtXwvYEUhizYzs6ZpcM8/Iq6JiF4RUUrmh9wnI+J8YCYwJuk2Bng0eT4TGC2pg6Q+QBkwu+CVm5lZ3ppySefJwAxJY4FlwCiAiFgoaQbwOrAFmOAjfczMWpdGhX9EPA08nTxfCwyvp98kYFITazMzsyLxGb5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZCDn8zsxRy+JuZpZDD38wshRz+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaVQg9fzl9QReBbokPR/MCKul3QQMB0oBZYC34uID5JlrgHGAluBKyLiT0Wp3iylxla+3Ohl7r5ocBEqsT1VLnv+m4FvRMSxwCBghKQTgYnArIgoA2Yl00jqT+Z2j0cDI4ApktoVoXYzM8tTLvfwjYjYkEzunfwFMBKoStqrgLOS5yOBaRGxOSKWAIuB4wtZtJmZNU1OY/6S2kmaD6wGnoiIl4BDImIlQPJ4cNK9J7A8a/HqpG3ndY6TNEfSnJqamia8BTMza6ycwj8itkbEIKAXcLykAbvprrpWUcc6KyKiPCLKS0pKcirWzMwKo1FH+0TEOjI3cB8BrJLUHSB5XJ10qwZ6Zy3WC1jR1ELNzKxwcjnapwT4LCLWSfoC8E3gl8BMYAwwOXl8NFlkJvCApFuBHkAZMLsIte+R8jlKw8ys0BoMf6A7UJUcsbMXMCMiHpP0AjBD0lhgGTAKICIWSpoBvA5sASZExNbilG9mZvloMPwj4lXgy3W0rwWG17PMJGBSk6szM7Oi8Bm+ZmYp5PA3M0shh7+ZWQo5/M3MUsjhb2aWQg5/M7MUcvibmaWQw9/MLIUc/mZmKeTwNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlUC538uoN3At8EfgcqIiI2yQdBEwHSoGlwPci4oNkmWuAscBW4IqI+FNRqjeznOV7F7m7Lxpc4EqsNchlz38L8LOI6AecCEyQ1B+YCMyKiDJgVjJNMm80cDSZe/1OSe4CZmZmrUSD4R8RKyNiXvL8I+ANoCcwEqhKulUBZyXPRwLTImJzRCwBFgPHF7huMzNrgkaN+UsqJXNLx5eAQyJiJWQ+IICDk249geVZi1UnbTuva5ykOZLm1NTU5FG6mZnlK+fwl7Q/8AfgJxGxfndd62iLXRoiKiKiPCLKS0pKci3DzMwKIKfwl7Q3meC/PyIeSppXSeqezO8OrE7aq4HeWYv3AlYUplwzMyuEBsNfkoC7gTci4tasWTOBMcnzMcCjWe2jJXWQ1AcoA2YXrmQzM2uqBg/1BE4CLgBekzQ/absWmAzMkDQWWAaMAoiIhZJmAK+TOVJoQkRsLXThZmaWvwbDPyKeo+5xfIDh9SwzCZjUhLrMzKyIfIavmVkKOfzNzFIolzF/q0e+p8ubmbU07/mbmaWQw9/MLIUc/mZmKeTwNzNLIYe/mVkKOfzNzFLI4W9mlkIOfzOzFHL4m5mlkMPfzCyFHP5mZink8DczS6Fc7uR1j6TVkhZktR0k6QlJbyePXbLmXSNpsaRFkk4tVuFmZpa/XPb8K4ERO7VNBGZFRBkwK5lGUn9gNHB0sswUSe0KVq2ZmRVEg+EfEc8C7+/UPBKoSp5XAWdltU+LiM0RsQRYDBxfmFLNzKxQ8r2e/yERsRIgIlZKOjhp7wm8mNWvOmkzsz1UPvetuPuiwUWoxAqp0D/41nWv36izozRO0hxJc2pqagpchpmZ7U6+4b9KUneA5HF10l4N9M7q1wtYUdcKIqIiIsojorykpCTPMszMLB/5hv9MYEzyfAzwaFb7aEkdJPUByoDZTSvRzMwKrcExf0lTgWFAN0nVwPXAZGCGpLHAMmAUQEQslDQDeB3YAkyIiK1Fqt3MzPLUYPhHxHn1zBpeT/9JwKSmFNUSfDN2M0sTn+FrZpZCDn8zsxRy+JuZpZDD38wshfI9w9fMrF4+K7j1856/mVkKOfzNzFLI4W9mlkIe8zezViHfEy39W0F+vOdvZpZCDn8zsxRqc8M+vkaPmVnDvOdvZpZCDn8zsxRy+JuZpZDD38wshYr2g6+kEcBtQDvgtxExuVivZWbp5esI5acoe/6S2gF3AqcB/YHzJPUvxmuZmVnjFWvY53hgcUS8GxGfAtOAkUV6LTMza6RiDfv0BJZnTVcDJ2R3kDQOGJdMbpa0oEi17Gm6AWtauohWwttiO2+L7Zq8Le65uECVtLwj812wWOGvOtpih4mICqACQNKciCgvUi17FG+L7bwttvO22M7bYjtJc/JdtljDPtVA76zpXsCKIr2WmZk1UrHC/2WgTFIfSfsAo4GZRXotMzNrpKIM+0TEFkmXAX8ic6jnPRGxcDeLVBSjjj2Ut8V23hbbeVts522xXd7bQhHRcC8zM2tTfIavmVkKOfzNzFKoWcNf0ghJiyQtljSxjvmSdHsy/1VJxzVnfc0ph23xT8k2eFXS85KObYk6m0ND2yKr32BJWyWd05z1NadctoWkYZLmS1oo6ZnmrrG55PD/SGdJ/yHpv5Nt0XaO3s8i6R5Jq+s7Fyrv3IyIZvkj88PvO8BhwD7AfwP9d+pzOvBfZM4TOBF4qbnqa86/HLfFEKBL8vy0NG+LrH5PAv8POKel627BfxcHAq8DX0qmD27pultwW1wL/DJ5XgK8D+zT0rUXYVsMBY4DFtQzP6/cbM49/1wu+TASuDcyXgQOlNS9GWtsLg1ui4h4PiI+SCZfJHOuRFuU66VALgf+AKxuzuKaWS7b4h+BhyJiGUBEtNXtkcu2CKCTJAH7kwn/Lc1bZvFFxLNk3lt98srN5gz/ui750DOPPm1BY9/nWDKf7G1Rg9tCUk/gfwB3NWNdLSGXfxd9gS6SnpY0V9KFzVZd88plW9wB9CNzAulrwI8j4vPmKa9VySs3m/Mevg1e8iHHPm1Bzu9T0ilkwv/kolbUcnLZFr8Gro6IrZmdvDYrl23RHvgKMBz4AvCCpBcj4q1iF9fMctkWpwLzgW8AhwNPSPpLRKwvcm2tTV652Zzhn8slH9JyWYic3qekY4DfAqdFxNpmqq255bItyoFpSfB3A06XtCUiHmmWCptPrv+PrImIj4GPJT0LHAu0tfDPZVtcDEyOzMD3YklLgKOA2c1TYquRV24257BPLpd8mAlcmPx6fSLwYUSsbMYam0uD20LSl4CHgAva4F5dtga3RUT0iYjSiCgFHgTGt8Hgh9z+H3kU+Jqk9pL2JXO13Deauc7mkMu2WEbmGxCSDiFzhct3m7XK1iGv3Gy2Pf+o55IPki5N5t9F5kiO04HFwEYyn+xtTo7b4jqgKzAl2ePdEm3wSoY5botUyGVbRMQbkv4IvAp8TuYueW3ucug5/rv4F6BS0mtkhj6ujog2d9lrSVOBYUA3SdXA9cDe0LTc9OUdzMxSyGf4mpmlkMPfzCyFHP5mZink8DczSyGHv5lZCjn8zcxSyOFvZpZC/x8X3zE2UQe2DwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# 收敛步数\n",
    "m = 1000\n",
    "# 迭代步数\n",
    "n = 10000\n",
    "\n",
    "# 建议分布\n",
    "proposal_dist = ProposalDistribution()\n",
    "# 接受分布\n",
    "accepted_dist = AcceptedDistribution()\n",
    "\n",
    "metropolis_hastings = MetropolisHastings(proposal_dist, accepted_dist, m, n)\n",
    "\n",
    "# 使用 Metropolis-Hastings 算法进行求解\n",
    "samples, dist_mean, dist_var = metropolis_hastings.solve()\n",
    "print(\"均值:\", dist_mean)\n",
    "print(\"方差:\", dist_var)\n",
    "\n",
    "# 对结果进行可视化\n",
    "metropolis_hastings.visualize(samples, bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题19.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;设某试验可能有五种结果，其出现的概率分别为\n",
    "$$\n",
    "\\frac{\\theta}{4}+\\frac{1}{8}, \\quad \\frac{\\theta}{4}, \\quad \\frac{\\eta}{4}, \\quad \\frac{\\eta}{4}+\\frac{3}{8}, \\quad \\frac{1}{2}(1-\\theta-\\eta)\n",
    "$$\n",
    "模型含有两个参数 $\\theta$ 和 $\\eta$，都介于 0 和 1 之间。现有 22 次试验结果的观测值为\n",
    "$$\n",
    "y = (y_1, y_2, y_3, y_4, y_5) = (14, 1, 1, 1, 5)\n",
    "$$\n",
    "其中 $y_i$ 表示 22 次试验中第 $i$ 个结果出现的次数，$i = 1,2,\\cdots,5$。试用吉布斯抽样估计参数 $\\theta$ 和 $\\eta$ 的均值和方差。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解题思路**\n",
    "\n",
    "1. 给出吉布斯抽样算法\n",
    "2. 定义目标概率的密度函数和各参数的满条件分布\n",
    "3. 自编程实现吉布斯抽样算法进行求解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：吉布斯抽样算法**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第19.5.2节的吉布斯抽样算法的描述：\n",
    "\n",
    "> **算法19.3(吉布斯抽样)**  \n",
    "> 输入：目标概率分布的密度函数 $p(x)$，函数 $f(x)$；  \n",
    "> 输出：$p(x)$ 的随机样本 $x_{m+1}, x_{m+2}, \\cdots, x_{n}$，函数样本均值 $f_{m n}$；  \n",
    "> 参数：收敛步数 $m$，迭代步数 $n$。  \n",
    "> （1）初始化。给出初始样本 $x^{(0)} = \\left( x_1^{(0)}, x_2^{(0)}, \\cdots, x_k^{(0)} \\right)^T$。  \n",
    "> （2）对 $i$ 循环执行  \n",
    "> 设第 $(i-1)$ 次迭代结束时的样本为 $x^{(i-1)} = \\left( x_1^{(i-1)}, x_2^{(i-1)}, \\cdots, x_k^{(i-1)} \\right)^T$，则第 $i$ 次迭代进行如下几步操作：\n",
    "> $$ \\left \\{ \n",
    "\\begin{array}{l}\n",
    "\\text {(1) 由满条件分布 } p \\left( x_1 | x_2^{(i-1)}, \\cdots, x_k^{(i-1)} \\right) \\text { 抽取 } x_1^{(i)} \\\\\n",
    "\\quad \\quad \\vdots \\\\\n",
    "\\text {(j) 由满条件分布 } p \\left( x_j | x_1^{(i)}, \\cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \\cdots, x_k^{(i-1)} \\right) \\text{ 抽取 } x_j^{(i)} \\\\\n",
    "\\quad \\quad \\vdots \\\\\n",
    "\\text {(k) 由满条件分布 } p \\left( x_k | x_1^{(i)}, \\cdots, x_{k-1}^{(i)} \\right) \\text { 抽取 } x_k^{(i)}\n",
    "\\end{array}\n",
    "\\right. $$\n",
    "> 得到第 $i$ 次迭代值 $x^{(i)} = \\left( x_1^{(i)}, x_2^{(i)}, \\cdots, x_k^{(i)} \\right)^T$。  \n",
    "> （3）得到样本集合\n",
    "> $$ \\{ x^{(m+1)}, x^{(m+2)}, \\cdots, x^{(n)} \\} $$\n",
    "> （4）计算\n",
    "> $$ f_{m n} = \\frac{1}{n-m} \\sum_{i=m+1}^n f (x^{(i)}) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：定义目标概率的密度函数和各参数的满条件分布**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据题意，可知目标概率的分布函数为\n",
    "$$\n",
    "\\begin{aligned}\n",
    "P\n",
    "&= (p_1, p_2, p_3, p_4, p_5) \\\\\n",
    "&= (\\frac{\\theta}{4}+\\frac{1}{8},\\ \\frac{\\theta}{4}, \\ \\frac{\\eta}{4}, \\ \\frac{\\eta}{4}+\\frac{3}{8}, \\ \\frac{1}{2}(1-\\theta-\\eta))\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;可得：\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\displaystyle \\frac{\\theta}{4}+\\frac{1}{8} > 0 \\\\\n",
    "\\displaystyle \\frac{\\theta}{4} > 0 \\\\\n",
    "\\displaystyle \\frac{\\eta}{4} > 0 \\\\\n",
    "\\displaystyle \\frac{\\eta}{4}+\\frac{3}{8} > 0 \\\\\n",
    "\\displaystyle \\frac{1}{2}(1-\\theta-\\eta) > 0 \\\\\n",
    "\\displaystyle \\frac{\\theta}{4}+\\frac{1}{8} + \\frac{\\theta}{4} + \\frac{\\eta}{4} + \\frac{\\eta}{4}+\\frac{3}{8} + \\frac{1}{2}(1-\\theta-\\eta) = 1 \\\\\n",
    "0 < \\theta < 1 \\\\\n",
    "0 < \\eta < 1 \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;求解得到：\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\theta + \\eta < 1 \\\\\n",
    "0 < \\theta < 1 \\\\\n",
    "0 < \\eta < 1 \\\\\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;则参数 $\\theta$ 和 $\\eta$ 的满分布条件概率为\n",
    "$$\n",
    "P (\\theta | y, \\eta) = \\left( \\left( \\frac{\\theta}{4}+\\frac{1}{8} \\right)^{14}, \\ \\frac{\\theta}{4}, \\ \\frac{\\eta}{4}, \\ \\frac{\\eta}{4}+\\frac{3}{8}, \\ \\left( \\frac{1}{2} (1-\\theta-\\eta) \\right)^5 \\right), \\quad \\theta \\in (0, 1-\\eta) \\\\\n",
    "P (\\eta | y, \\theta) = \\left( \\left( \\frac{\\theta}{4}+\\frac{1}{8} \\right)^{14}, \\ \\frac{\\theta}{4}, \\ \\frac{\\eta}{4}, \\ \\frac{\\eta}{4}+\\frac{3}{8}, \\ \\left( \\frac{1}{2} (1-\\theta-\\eta) \\right)^5 \\right), \\quad \\eta \\in (0, 1-\\theta)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：编程实现吉布斯抽样算法进行求解**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "class GibbsSampling:\n",
    "    def __init__(self, target_dist, j, m=1e4, n=1e5):\n",
    "        \"\"\"\n",
    "        Gibbs Sampling 算法\n",
    "\n",
    "        :param target_dist: 目标分布\n",
    "        :param j: 变量维度\n",
    "        :param m: 收敛步数\n",
    "        :param n: 迭代步数\n",
    "        \"\"\"\n",
    "        self.target_dist = target_dist\n",
    "        self.j = j\n",
    "        self.m = int(m)\n",
    "        self.n = int(n)\n",
    "\n",
    "    def solve(self):\n",
    "        \"\"\"\n",
    "        Gibbs Sampling 算法求解\n",
    "        \"\"\"\n",
    "        # (1) 初始化\n",
    "        all_samples = np.zeros((self.n, self.j))\n",
    "        # 任意选择一个初始值\n",
    "        x_0 = np.random.random(self.j)\n",
    "        # x_0 = np.array([0.5, 0.5])\n",
    "        # (2) 循环执行\n",
    "        for i in range(self.n):\n",
    "            x = x_0 if i == 0 else all_samples[i - 1]\n",
    "            # 满条件分布抽取\n",
    "            for k in range(self.j):\n",
    "                x[k] = self.target_dist.sample(x, k)\n",
    "            all_samples[i] = x\n",
    "        # (3) 得到样本集合\n",
    "        samples = all_samples[self.m:]\n",
    "        # (4) 计算函数样本均值\n",
    "        dist_mean = samples.mean(0)\n",
    "        dist_var = samples.var(0)\n",
    "        return samples[self.m:], dist_mean, dist_var\n",
    "\n",
    "    @staticmethod\n",
    "    def visualize(samples, bins=50):\n",
    "        \"\"\"\n",
    "        可视化展示\n",
    "        :param samples: 抽取的随机样本集合\n",
    "        :param bins: 频率直方图的分组个数\n",
    "        \"\"\"\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.set_title('Gibbs Sampling')\n",
    "        ax.hist(samples[:, 0], bins, alpha=0.7, label='$\\\\theta$')\n",
    "        ax.hist(samples[:, 1], bins, alpha=0.7, label='$\\\\eta$')\n",
    "        ax.set_xlim(0, 1)\n",
    "        ax.legend()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TargetDistribution:\n",
    "    \"\"\"\n",
    "    目标概率分布\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self):\n",
    "        # 联合概率值过小，可对建议分布进行放缩\n",
    "        self.c = self.__select_prob_scaler()\n",
    "\n",
    "    def sample(self, x, k=0):\n",
    "        \"\"\"\n",
    "        使用接受-拒绝方法从满条件分布中抽取新的分量 x_k\n",
    "        \"\"\"\n",
    "        theta, eta = x\n",
    "        if k == 0:\n",
    "            while True:\n",
    "                new_theta = np.random.uniform(0, 1 - eta)\n",
    "                alpha = np.random.uniform()\n",
    "                if (alpha * self.c) < self.__prob([new_theta, eta]):\n",
    "                    return new_theta\n",
    "        elif k == 1:\n",
    "            while True:\n",
    "                new_eta = np.random.uniform(0, 1 - theta)\n",
    "                alpha = np.random.uniform()\n",
    "                if (alpha * self.c) < self.__prob([theta, new_eta]):\n",
    "                    return new_eta\n",
    "\n",
    "    def __select_prob_scaler(self):\n",
    "        \"\"\"\n",
    "        选择合适的建议分布放缩尺度\n",
    "        \"\"\"\n",
    "        prob_list = []\n",
    "        step = 1e-3\n",
    "        for theta in np.arange(step, 1, step):\n",
    "            for eta in np.arange(step, 1 - theta + step, step):\n",
    "                prob = self.__prob((theta, eta))\n",
    "                prob_list.append(prob)\n",
    "        searched_max_prob = max(prob_list)\n",
    "        upper_bound_prob = searched_max_prob * 10\n",
    "        return upper_bound_prob\n",
    "\n",
    "    @staticmethod\n",
    "    def __prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        theta = x[0]\n",
    "        eta = x[1]\n",
    "        p1 = (theta / 4 + 1 / 8) ** 14\n",
    "        p2 = theta / 4\n",
    "        p3 = eta / 4\n",
    "        p4 = (eta / 4 + 3 / 8)\n",
    "        p5 = 1 / 2 * (1 - theta - eta) ** 5\n",
    "        p = (p1 * p2 * p3 * p4 * p5)\n",
    "        return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta均值：0.5201289483641778, theta方差：0.017691128644918685\n",
      "eta均值：0.12203712909117162, eta方差：0.0064478491785926505\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYoUlEQVR4nO3de7SddX3n8feHcAkIlECAIoEmYLwA3gOKtg4drCBooWvJELAKli6WFS8zaqvY6bi0Zqm1M16mRZuKAqsoMlwGlNHhYpFx5BZqBAKiQZScIZIQQClKJPidP/bDk93jSc7J2fvss3Pyfq111nme3/N7nue3fznZn/17bjtVhSRJANtNdwMkScPDUJAktQwFSVLLUJAktQwFSVLLUJAktQwFbRWSfC7JXzXTRyUZ2Uzd85J8ZHCtm1pJTk/y7a75f01y0HS2STOXoaChkGRxkpuTPJ5kTTP9tiQBqKq3VtVfD7hNhya5OskjSR5NcluS4wbZhrFU1a5V9aPpbodmJkNB0y7Je4BPA58AfhvYF3gr8Epgx2ls2leBa5r27AO8E/j5NLZHmnKGgqZVkt8CPgy8raouqarHquO7VfXGqlrf1PuNQ0JJPpDkoSQ/TvLGUZuem+SaJI8l+VaS32nWSZJPNqORnyW5PclhY7RrLrAA+Meq+lXz83+r6tvN8jlJvpZkbTOS+FqSeV3rX5/kI0m+0xzu+WqSvZJcmOTnSW5NMr+rfiV5Z5IfNa/pE0nG/P/Z1H1WV7/8fZKrmtd6c5KDu+q+Jsk9zWs9p+mLP92SfyNtWwwFTbcjgZ2AK7Zwvd8G5gL7A6cBS5M8p2v5G4G/buosBy5syl8DvAp4NrAHcDKwboztrwNWAv+U5MQk+45avh3wReB3gAOBXwJ/N6rOYuBNTRsPBm5s1tkTuBv44Kj6fwQsAl4CnAD8yaZf/r9xCvAhYE7T5iXQBtslwNnAXsA9wCsmuE1towwFTbe5wENVteHpgubT9aNJfpnkVZtZ96+qan1VfQu4CvgPXcuuqqobmpHGXwJHJjkAeBLYDXgukKq6u6pWj95wdR4K9vvAj4H/CqxOckOShc3ydVV1aVX9oqoeo/NG/O9GbeaLVXVvVf0M+Dpwb1Vd27zW/wG8eFT9j1fVw1V1P/ApOm/2E3FZVd3SbPdC4EVN+XHAiqq6rFn2GeCnE9ymtlGGgqbbOjqHerZ/uqCqXlFVezTLNvU3+khVPd41/xPgmV3zq7q296/Aw8Azq+qbdD7R/z3wYJKlSXYfawdVNVJVb6+qg+mMCB4HLgBIskuSf0jykyQ/B24A9kgyq2sTD3ZN/3KM+V1H7XJV1/To17M53W/0v+ja7jP5t/1QwCav2pLAUND0uxFYT+dwyZaYk+QZXfMHAg90zR/w9ESSXekcsnkAoKo+U1UvBQ6lcxjpz8fbWVWtohMkT59/eA/wHOBlVbU7nUNSANnC19HtgK7p0a9nMlYD3ec50j0vjcVQ0LSqqkfpHA8/J8kbkuyaZLskLwKesdmV4UNJdkzye8Dr6BySedpxSX43yY50zi3cXFWrkhye5GVJdqDzyf8J4KnRG25OJH8oybOa9sylc4z/pqbKbnQ+7T+aZE9+8/zAZPx5s98DgHcBX+lxe1cBz2/OiWwPnEXnXIy0SYaCpl1V/Q3wbuAvgDV0DrP8A/A+4DubWO2nwCN0Pk1fCLy1qr7ftfxLdN6oHwZeSufEM8DuwD826/6EziGqvx1j+78C5gPX0rkM9U46I5rTm+WfAnYGHqITFN+Y6OvdjCuA2+icGL8KOLeXjVXVQ8BJwN/QeZ2HAMvovA5pTPFLdqTpl6SAhVW1cgr3sR2dcwpvrKp/nqr9aOvmSEGawZIck2SPJDsBH6BzzuOmcVbTNsxQkGa2I4F76Rzmej1wYlX9cnqbpGHm4SNJUsuRgiSptf34VabX3Llza/78+dPdDEnaqtx2220PVdXeW7re0IfC/PnzWbZs2XQ3Q5K2Kkl+Mpn1PHwkSWoZCpKklqEgSWoN/TkFSRq0J598kpGREZ544onpbsq4Zs+ezbx589hhhx36sj1DQZJGGRkZYbfddmP+/Pk0XxM+lKqKdevWMTIywoIFC/qyTQ8fSdIoTzzxBHvttddQBwJAEvbaa6++jmgMBUkaw7AHwtP63U5DQZLUMhQkSa1t40Tzl04ev86pvX7JlSRt/baNUJCkHpxx3q193d65px8+oXpPPfUU7373u7n22mvZbrvtuOKKKzjooIP62pbRPHwkSUPqox/9KAcddBArVqzgne98J+ecc86U79ORgiQNoccff5zLL7+c2267DYAFCxZw1VVXTfl+DQVJGkLXXnstq1at4kUvehEADz/8MK9+9aunfL/jHj5K8oUka5Lc2VX2iSTfT3J7ksuT7NG17OwkK5Pck+SYrvKXJrmjWfaZbC0XAUvSNFi+fDkf/vCHWb58OcuXL+c1r3lNGxBTaSLnFM4Djh1Vdg1wWFW9APgBcDZAkkOAxcChzTrnJJnVrPNZ4ExgYfMzepuSpMYjjzzCLrvsAsCGDRu4+uqref3rXz/l+x03FKrqBuDhUWVXV9WGZvYmYF4zfQJwUVWtr6r7gJXAEUn2A3avqhur86XQFwAn9uk1SNKM8+xnP5ubbroJgE9+8pMcf/zxfXu+0eb045zCnwBPX+S/P52QeNpIU/ZkMz26fExJzqQzquDAAw/sQxMlafImeglpP51yyim89rWv5VnPehZHHnkkS5cuHch+ewqFJH8JbAAufLpojGq1mfIxVdVSYCnAokWLNllPkmaqOXPmtCOFQZp0KCQ5DXgdcHRzSAg6I4ADuqrNAx5oyueNUS5JGiKTunktybHA+4A/rKpfdC26ElicZKckC+icUL6lqlYDjyV5eXPV0ZuBK3psuySpz8YdKST5MnAUMDfJCPBBOlcb7QRc01xZelNVvbWqViS5GLiLzmGls6rqqWZTf0bnSqadga83P8PD5yNJ0vihUFWnjFF87mbqLwGWjFG+DDhsi1onSRoon30kSWoZCpKkls8+kqTxTOSc45YY4vOTjhQkSS1HCpI0pD7+8Y9z7733smbNGr73ve9x1lln8d73vndK9+lIQZKG1B133MH69eu59NJLufrqq7ngggumfJ+OFCRpSN1+++1cfvnlzJo1i1mzZrHnnntO+T4dKUjSEHryySd56KGHOPjgg4FOQDz/+c+f8v0aCpI0hO655x6e97zntfPLly/nhS984ZTv18NHkjSeabiE9I477uAFL3hBO798+XKOP/74Kd+voSBJQ+iUU07hlFM2PmXozjvv5LDDpv5JQTMjFPp9Y4kkDZFHH32UHXfckZ133nnK9+U5BUkacnvssQd33XXXQPZlKEiSWoaCJKllKEjSGDZ+y/Bw63c7DQVJGmX27NmsW7du6IOhqli3bh2zZ8/u2zZnxtVHktRH8+bNY2RkhLVr1053U8Y1e/Zs5s2b17ftGQqSNMoOO+zAggULprsZ08LDR5KklqEgSWoZCpKklqEgSWoZCpKk1rihkOQLSdYkubOrbM8k1yT5YfN7Tteys5OsTHJPkmO6yl+a5I5m2WeSpP8vR5LUi4mMFM4Djh1V9n7guqpaCFzXzJPkEGAxcGizzjlJZjXrfBY4E1jY/IzepiRpmo17n0JV3ZBk/qjiE4CjmunzgeuB9zXlF1XVeuC+JCuBI5L8GNi9qm4ESHIBcCLw9Z5fgTSNzjjv1kmve+7ph/exJVJ/TPacwr5VtRqg+b1PU74/sKqr3khTtn8zPbpckjRE+n2ieazzBLWZ8rE3kpyZZFmSZVvDbeaSNFNM9jEXDybZr6pWJ9kPWNOUjwAHdNWbBzzQlM8bo3xMVbUUWAqwaNGi4Xki1US+4W0avstVkvplsiOFK4HTmunTgCu6yhcn2SnJAjonlG9pDjE9luTlzVVHb+5aR5I0JMYdKST5Mp2TynOTjAAfBD4GXJzkDOB+4CSAqlqR5GLgLmADcFZVPdVs6s/oXMm0M50TzJ5klqQhM5Grj07ZxKKjN1F/CbBkjPJlwGFb1DppBpvslUtetaSp5B3NkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqSWoSBJahkKkqTWZB+IJ80ovXwvgjSTOFKQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy8dcaMbwURVS7wwFaSvTS/ide/rhfWyJZiIPH0mSWoaCJKnVUygk+U9JViS5M8mXk8xOsmeSa5L8sPk9p6v+2UlWJrknyTG9N1+S1E+TDoUk+wPvBBZV1WHALGAx8H7guqpaCFzXzJPkkGb5ocCxwDlJZvXWfElSP/V6onl7YOckTwK7AA8AZwNHNcvPB64H3gecAFxUVeuB+5KsBI4AbuyxDcPlSyePX+fUr0x9OyRpEiY9Uqiq/wf8LXA/sBr4WVVdDexbVaubOquBfZpV9gdWdW1ipCmTJA2JXg4fzaHz6X8B8EzgGUn+eHOrjFFWm9j2mUmWJVm2du3ayTZRkrSFejnR/GrgvqpaW1VPApcBrwAeTLIfQPN7TVN/BDiga/15dA43/YaqWlpVi6pq0d57791DEyVJW6KXULgfeHmSXZIEOBq4G7gSOK2pcxpwRTN9JbA4yU5JFgALgVt62L8kqc8mfaK5qm5OcgnwL8AG4LvAUmBX4OIkZ9AJjpOa+iuSXAzc1dQ/q6qe6rH9kqQ+6unqo6r6IPDBUcXr6Ywaxqq/BFjSyz4lSVPHO5olSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSa1ev3lN0lbkjPNundR6555+eJ9bomHlSEGS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1PI+henwpZPHr3PqV6a+HZI0iiMFSVLLUJAktTx8pKEz2UcxSOqdIwVJUqunUEiyR5JLknw/yd1JjkyyZ5Jrkvyw+T2nq/7ZSVYmuSfJMb03X5LUT72OFD4NfKOqngu8ELgbeD9wXVUtBK5r5klyCLAYOBQ4Fjgnyawe9y9J6qNJh0KS3YFXAecCVNWvqupR4ATg/Kba+cCJzfQJwEVVtb6q7gNWAkdMdv+SpP7rZaRwELAW+GKS7yb5fJJnAPtW1WqA5vc+Tf39gVVd6480ZZKkIdFLKGwPvAT4bFW9GHic5lDRJmSMshqzYnJmkmVJlq1du7aHJkqStkQvoTACjFTVzc38JXRC4sEk+wE0v9d01T+ga/15wANjbbiqllbVoqpatPfee/fQREnSlph0KFTVT4FVSZ7TFB0N3AVcCZzWlJ0GXNFMXwksTrJTkgXAQuCWye5fktR/vd689g7gwiQ7Aj8C3kInaC5OcgZwP3ASQFWtSHIxneDYAJxVVU/1uH9JUh/1FApVtRxYNMaiozdRfwmwpJd9SpKmjnc0S5JahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJavT4QT9I24Izzbp30uueefngfW6Kp5khBktQa/pHCwz+CL5083a2QpG2CIwVJUstQkCS1hv/w0bZqvENmp35lMO2QtE1xpCBJahkKkqSWoSBJanlOQVOil5udJE0fRwqSpJahIElqGQqSpJahIElq9RwKSWYl+W6SrzXzeya5JskPm99zuuqenWRlknuSHNPrviVJ/dWPkcK7gLu75t8PXFdVC4HrmnmSHAIsBg4FjgXOSTKrD/uXJPVJT6GQZB5wPPD5ruITgPOb6fOBE7vKL6qq9VV1H7ASOKKX/UuS+qvXkcKngL8Aft1Vtm9VrQZofu/TlO8PrOqqN9KU/YYkZyZZlmTZ2sfW99hESdJETToUkrwOWFNVt010lTHKaqyKVbW0qhZV1aK9d9tpsk2UJG2hXu5ofiXwh0mOA2YDuyf5J+DBJPtV1eok+wFrmvojwAFd688DHuhh/5KkPpv0SKGqzq6qeVU1n84J5G9W1R8DVwKnNdVOA65opq8EFifZKckCYCFwy6RbLknqu6l49tHHgIuTnAHcD5wEUFUrklwM3AVsAM6qqqemYP+SpEnqSyhU1fXA9c30OuDoTdRbAizpxz4lSf3nHc2SpJahIElqGQqSpJahIElqGQqSpJZfx7m1+tLJ49c59StT3w5JM4ojBUlSy1CQJLUMBUlSy1CQJLUMBUlSy6uPtFlnnHfrdDdB0gA5UpAktRwpSJpSkx1tnnv64X1uiSbCkYIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJaXpI6k/l4bUlbyJGCJKllKEiSWoaCJKk16VBIckCSf05yd5IVSd7VlO+Z5JokP2x+z+la5+wkK5Pck+SYfrwASVL/9DJS2AC8p6qeB7wcOCvJIcD7geuqaiFwXTNPs2wxcChwLHBOklm9NF6S1F+TDoWqWl1V/9JMPwbcDewPnACc31Q7HzixmT4BuKiq1lfVfcBK4IjJ7l+S1H99OaeQZD7wYuBmYN+qWg2d4AD2aartD6zqWm2kKRtre2cmWZZk2drH1vejiZKkCeg5FJLsClwK/Meq+vnmqo5RVmNVrKqlVbWoqhbtvdtOvTZRkjRBPYVCkh3oBMKFVXVZU/xgkv2a5fsBa5ryEeCArtXnAQ/0sn9JUn/1cvVRgHOBu6vqv3UtuhI4rZk+Dbiiq3xxkp2SLAAWArdMdv+SpP7r5TEXrwTeBNyRZHlT9gHgY8DFSc4A7gdOAqiqFUkuBu6ic+XSWVX1VA/7lyT12aRDoaq+zdjnCQCO3sQ6S4Alk92nJmf5qkc3vfDjm79d5B3Af9/3I31tj6Th5R3NkqSWoSBJahkKkqSWoSBJahkKkqSW37wmaSidcd6tk1733NMP72NLti2OFCRJLUNBktQyFCRJLc8pbEU2e2fyFHrHg/953Dre9SzNDI4UJEktQ0GS1PLwkfpivENMHl6Stg6OFCRJLUNBktQyFCRJLUNBktTyRPOATde9BpI0EY4UJEktQ0GS1PLwkQbCR2VokCb72G0fue1IQZLUxZGChoajCWn6GQqT5FVEkmaigYdCkmOBTwOzgM9X1ccG3QZtvRxNSFMrVTW4nSWzgB8AfwCMALcCp1TVXZtaZ9FBe9ayj/zBlLTHT/vbLoND/TZsJ6mT3FZVi7Z0vUGPFI4AVlbVjwCSXAScAGwyFKSpMJERRz8YPtraDHqk8Abg2Kr602b+TcDLqurto+qdCZzZzB4G3DmwRg63ucBD092IIWFfbGRfbGRfbPScqtptS1ca9EghY5T9RipV1VJgKUCSZZMZAs1E9sVG9sVG9sVG9sVGSZZNZr1B36cwAhzQNT8PeGDAbZAkbcKgQ+FWYGGSBUl2BBYDVw64DZKkTRjo4aOq2pDk7cD/pnNJ6heqasU4qy2d+pZtNeyLjeyLjeyLjeyLjSbVFwM90SxJGm4++0iS1DIUJEmtoQiFJMcmuSfJyiTvH2N5knymWX57kpdMRzsHYQJ98camD25P8p0kL5yOdg7CeH3RVe/wJE8198HMSBPpiyRHJVmeZEWSbw26jYMygf8jv5Xkq0m+1/TFW6ajnYOQ5AtJ1iQZ816uSb13VtW0/tA54XwvcBCwI/A94JBRdY4Dvk7nPoeXAzdPd7unsS9eAcxppl+7LfdFV71vAv8LeMN0t3sa/y72oPNkgAOb+X2mu93T2BcfAD7eTO8NPAzsON1tn6L+eBXwEuDOTSzf4vfOYRgptI++qKpfAU8/+qLbCcAF1XETsEeS/Qbd0AEYty+q6jtV9UgzexOdez1moon8XQC8A7gUWDPIxg3YRPriVOCyqrofoKpman9MpC8K2C1JgF3phMKGwTZzMKrqBjqvb1O2+L1zGEJhf2BV1/xIU7aldWaCLX2dZ9D5FDATjdsXSfYH/gj43ADbNR0m8nfxbGBOkuuT3JbkzQNr3WBNpC/+DngenRtj7wDeVVW/Hkzzhs4Wv3cOw/cpTOTRFxN6PMYMMOHXmeT36YTC705pi6bPRPriU8D7quqpzofCGWsifbE98FLgaGBn4MYkN1XVD6a6cQM2kb44BlgO/HvgYOCaJP+nqn4+xW0bRlv83jkMoTCRR19sK4/HmNDrTPIC4PPAa6tq3YDaNmgT6YtFwEVNIMwFjkuyoar+50BaODgT/T/yUFU9Djye5AbghXQeVT+TTKQv3gJ8rDoH1VcmuQ94LnDLYJo4VLb4vXMYDh9N5NEXVwJvbs6kvxz4WVWtHnRDB2DcvkhyIHAZ8KYZ+Cmw27h9UVULqmp+Vc0HLgHeNgMDASb2f+QK4PeSbJ9kF+BlwN0DbucgTKQv7qczYiLJvsBzgB8NtJXDY4vfO6d9pFCbePRFkrc2yz9H58qS44CVwC/ofBKYcSbYF/8F2As4p/mEvKFm4FMhJ9gX24SJ9EVV3Z3kG8DtwK/pfKvhjHvk/AT/Lv4aOC/JHXQOn7yvqmbk47STfBk4CpibZAT4ILADTP6908dcSJJaw3D4SJI0JAwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktf4/JPk8AlOnbYcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 收敛步数\n",
    "m = 1e3\n",
    "# 迭代步数\n",
    "n = 1e4\n",
    "\n",
    "# 目标分布\n",
    "target_dist = TargetDistribution()\n",
    "\n",
    "# 使用 Gibbs Sampling 算法进行求解\n",
    "gibbs_sampling = GibbsSampling(target_dist, 2, m, n)\n",
    "\n",
    "samples, dist_mean, dist_var = gibbs_sampling.solve()\n",
    "\n",
    "print(f'theta均值：{dist_mean[0]}, theta方差：{dist_var[0]}')\n",
    "print(f'eta均值：{dist_mean[1]}, eta方差：{dist_var[1]}')\n",
    "\n",
    "# 对结果进行可视化\n",
    "GibbsSampling.visualize(samples, bins=20)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "273.188px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "fe245b48eec612fbee62cf5f1e4a377bd40a15fc2ce06ad16b3622723b03fe09"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
